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ABSTRACT 

We analyse the evaporation and condensation of spherical and cylindrical HI clouds of the cold 
neutral medium surrounded by the warm neutral medium. Because the interstellar medium 
including those two phases is well described as a thermally bistable fluid, it is useful to apply 
pattern formation theories to the dynamics of the interface between the two phases. Assuming 
isobaric evolution of fluids and a simple polynomial form of the heat-loss function, we show 
the curvature effects of the interface. We And that approximate solutions for spherical clouds 
are in good agreement with numerically obtained solutions. We extend our analysis to general 
curved interfaces taking into account the curvature effects explicitly. We And that the curvature 
effects always stabilise curved interfaces under assumptions such as isobaric evolution we 
adopt in this Letter. 

Key words: hydrodynamics - ISM: clouds - ISM: kinematics and dynamics - methods: 
analytical 


1 INTRODUCTION 


It is widely known that the interstellar medium (ISM) is well 
described as a thermally bistable fluid owing to radiative cool- 
ing and heating due to external radiation fields and c osmic rays 
^Field, Goldsmith & Habin^ll969l: l^fire et aljElol . Two sta¬ 
ble phases are called the warm neutral medium (WNM) with tem¬ 
perature T ~ 10'* K and the cold neutral medium (CNM) with 
T ~ 10*“^ K, respectively. Gas in an unstable phase with temper¬ 
ature between the WNM and CNM i s decompos e d into the tw o 
stable phases via thermal instability jpieljllQ^ iBalbuslIlQSfl) . 
These stable phases can be connected through interfaces in pressure 
equilibrium. It is important for unders tanding the behaviour of the 
ISM such as the interstellar turbulence iKovama & Inutsukal2002aL 


l20Q4t[&itsuk & Normanll2002J h^jAudi^^^ennebe^ l200^ ~and 


the evolution of galaxies (e.g. lMcKe^^^strikeil^^ ) to clarify 

the dynamics of the interface or front, which corresponds to the 
evap oration and condensation of low te mperature HI clouds . 

IZekdovich & PikeFneJ ^9^) and IPenston & BrownI jl97(]|) 
considered isobaric and steady phase-change in the bistable 
fluid assuming plane-parallel geometry. In this case, because 
the mass flux across the interface conserves, the velocity of 
fluids can be estimated as an eigenvalue of the energy equation. 
Then they pointed out the relationship between the motion 
of fronts, which determines the rate of evaporation or con- 
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densation of clouds, and external pressure, and the existence 
of saturation pressure at which a static front can exist. After 
that, many authors have investi g ated t he thermally bistable 
flow (e.g^_|Fem^^^hchekino^_^93|_|lfanebelle&Pera^ 

Ei pfikkiRege^^^plegeir^h^ir and lElnhlc^^egev&'shavM 

<1991 hereafter ERS92) treated with fluid equations in a more 
sophisticated manner in plane-parallel geometry. In the latter work, 
they formulated equations in the Lagrangian coordinate. Combined 
with the isobaric assumption, they derived a second-order ordinary 
differential equation in a steady state from the energy equation. 
Moreover, they discussed the behaviour of steady solutions from a 
pattern-theoretical point of view. 

It is challenging to extend those analyses to higher dimen- 
sions while it is much us eful for applying to realistic situations. 
iGraham & Langeil <1973h numerically computed isobaric flows 
in three dimensional spherical geometry for the first time. They 
pointed out the existence of a critical size of cloud s to avoid 
evaporation. Based on ERS92, IShaviv & RegevI <1994) argued the 
case of higher dimensions. However, it is difficult to extend the 
Lagrangian formulation in plane-parallel geometry provided by 
ERS92 to higher dimensional geometry in a straightforward man¬ 
ner. Therefore assuming a model equation similar to the Ginzburg- 
Landau (GL) equation, they derived the speed of frontal motion. 
Surprisingly, they have showed that the dependence of the frontal 
speed on time, or radii of clouds, is dependent on the dimen¬ 
sion of geometry. Besides they have claimed that their conclu- 
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sion is supported by numerical simu latio ns. It should be n oted that 
lAranson. Meerson & SasorovI fl995^ an eerso also dis¬ 

cussed the front curvature effects in a confined plasma, in which the 
boundary conditions are different from ours. Thereby their methods 
are inapplicable to the ISM in a straightforward way. 

Thus, the purpose of this paper is to reanalyse the frontal mo¬ 
tion in d-dimensional spheric ally symmet ric geometry, aided by 
pattern formation theories fe.g. lBravll994l . This paper is outlined 
as follows. In Section 2 we briefly review ERS92 to prepare the 
analysis in higher dimensional geometry. In Section 3 we show a 
systematic procedure to derive the curvature effects. In Section 4 
we discuss the dynamics of general curved fronts. In Section 5 we 
provide conclusions and discussion. 


2 LAGRANGIAN DESCRIPTION IN PLANE-PARALLEL 
GEOMETRY 


In the following, we assume isobaric evolution in which pressure 
is uniform over the whole system because the fluid motion is much 
slower than the sound speed. The basic fluid equations under the 
assumption of isobaric evolution are written as 

| + V.,v = 0 , ( 1 ) 

-v TL dr 

—^ -p^ + pr - V • k{t)vt = 0, (2) 

7 — 1 /i df 
7^ 

P = p-r, (3) 

p 

where the equation of motion is omitted because of the isobaric 
assumption, k(T’) is the conductivity, TZ the gas constant, p the 
mean molecular weight, and 7 the adiabatic index. The heat-loss 
function, £, is defined as pA — T, where A and T are the cooling 
and heating rates, respectively. Throughout this paper, we assume 
a single power-law of the conductivity, n{T) = noiT/To)°‘, and 
a = II 2 for neutral gas. The energy equation ^ under the iso¬ 
baric condition is derived by applying rds = d/i, where dp/p 
is neglected for the constant pressure, an d s and h are specific en¬ 
tro py and enthalpy, respective ly [see, e.g.. lPenston &BrowiJ<197(i) 
and lGraham & Lan gen J 1973^ 1. Dividing variables by characteristic 
values, we obtain dimensionless quantities, T = T/To, p = p/po, 
and p = p/po, and introducing characteristic length If and time- 
scale to, we obtain v = v/{I f/ to). Here we take to as the cooling 
time-scale. 


7 TZTo 
7 - 1 ppoAo ’ 


(4) 


where Aq is the cooling r ate at the ch aracteristic temperature, To, 
and IF as the Field length iFielJl96.5ll . 


It is useful to transform these equations into the Lagrangian 
description particularly in one dimensional case as shown in 
ERS92. Using a relation dm = pda;, where x is the dimension¬ 
less Eulerian coordinate normalised by If and m is the Lagrangian 
coordinate, the energy equation 0 is converted as follows, 

drf =-C+ pdmf°‘-^dmf, (9) 


where (r, m) are the independent dimensionless variables and dt 
and dm denote (d/dt) and {d/dm), respectively. Further transfor¬ 
mation introducing Z = T°‘ makes the above equation simpler, 

drZ^F{z)+pzPdlz\ ( 10 ) 

where = 1 — 1/a and F{Z) = —aZ^C. Note that this simple 
form is obtained only in the one dimensional case. ERS92 assumed 
a simple form of F{Z), which mimics the real heat-loss function, 

F{Z) = zP[dd{Z - Z 2 ) -{Z- Z 2 f - Tlogp], (11) 

where D is a positive constant. 

To obtain a travelling wave solution, we define x = ^ ~ 
CmT, where Cm is a constant which is a velocity on the Lagrangian 
coordinate. Then the partial differential equation llOt becomes an 
ordinary differential equation, 

pZ^Z” + CmZ' + F{Z) = 0, (12) 


where the prime means differentiation with respect to x- 

An analytic kink solution connecting the two stable phases is 
obtained when = 0 as shown in ERS92. The solution Zo is 


Zo{m) = Z 2 -\- Atanh 




(13) 


Unfortunate ly, _ as mentioned above and as shown in 

IShaviv & RegevI ^1994^ . it is difficult to extend it to higher dimen¬ 
sional cases. This can be easily seen if we take a d-dimensional 
(d 2) continuity, dm oc r'^~^Ax. Because this dependence 
on d produces a term explicitly depending on m, the second 
derivative dm cannot be regarded as V^. To avoid this diffi- 


^94 ) assum ed that a time-dependent GL 
Bray Il994h in the Lagrangian space is a 


culty, IShaviv & RegevI 
(TDGL) equation (e.g. 
good approximation for the evolution of a spherical cloud. To be 
free from such an assumption, it is fruitful to be back in the Eule¬ 
rian space to treat the fluid equations in higher dimensions. 


3 HIGHER DIMENSIONAL CASES: SPHERICAL AND 
CYLINDRICAL GEOMETRY 

In the Eulerian space (r, r), introducing X — obtain the 

continuity and energy equations. 


If = 


koTo 


(5) 


Using these quantities, the basic equations become the following 
dimensionless form, 

+ V • pv = 0. 

dr - 1 - - 

V-f T --;V • T“vr = 0, 

dr p 

p = pf, 


^-(l-ba)XV-v = 0, 
dr 


dX 1 + a . 


^ -b -xv^x. 

dr a P 


(14) 

(15) 


Transforming the coordinate to x = r — Rd{T) for a d-dimensional 

.d -1 „ 


( 6 ) 

spherical cloud. 


- RdX' = 

-vX' + {l + a)Xv' -b (1 -ba; 

vX' + 

(7) 

-RdX' = 

( 8 ) 


a 



+ ^x"x+ ^'^~^x'x. 


p r 


where r = t/to and C = £/poAo. 


(17) 
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where Rd — dRd/dr and F is given by ea. il 1> . The radius of the 
cloud, Rdir), is hereafter defined to be the position at X = X 2 = 
^(i+“)/“ ■pjjg velocity v is defined in the rest frame of the 
centre of the cloud. By defining Ud = v — Rd, which is the fluid 
velocity in the rest frame of the front, the above equations become 

UdX = (1 + q)Xm^ + (1 + a?)- {ud-\-Rd)X, (18) 

r 

UdX' = + -X''X + 

a p p r 

Because we take u = 0 at the cloud centre, Ud must be —Rd in 
the vicinity of the centre. In the plane-parallel case d = 1, using 
the relationship f = and pUd = —Cm, eQ.<17t 

can be reduced to eq.0, and therefore the kink solution ea. ilQt is 
easily proved to satisfy these equations when Ri = 0 . 

The above equations are non-linear, so it is hard to find an 
analytic solution. Therefore we simply assume what follows. One is 
that the structure of a solution is very similar to a one-dimensional 
(d = 1) solution. The second assumption is that the first derivative 
of the solution is sharply peaked at the interface, in other words, 
the width of the interface is very narrow compared to the scale of 
the cloud size. This corresponds to X' = 0 at r 7 ^ Rd- Under 
these assumptions, it is reasonable to substitute the first and second 
terms of the right-hand side of ea. <19> into uiX'. Thus we obtain 
an approximate form, 

Ud{R)=udR) + \'^^X2, ( 20 ) 

P 

where Ud{R) is the fluid velocity passing the front at r = in 
the front rest frame. X 2 in the last term emerges because we de¬ 
fine the position of the front at which X = X 2 . Thi s is similar 
to an equation discussed by lOraham & LangeJ il973h . Note that 
the above approximation corresponds to taking the first order in the 
expansion of the curvature term, 1/r = 1/R + 0{R~^), where 
1/77 is the ratio of the Field length If to the radius of the cloud 
^Aranson, Meerson & Sasorov|ll995ll . Thus this is valid only for 
7?> 1. 

FigQshows the fluid velocity at the front, Ud{R), against Rd 
for d = 2 and 3. The curves and symbols represent those given 
by the approximate solutions and numerical ones, respectively, for 
different values of pressure. Clearly the above approximate solution 
well agrees with the numerical solutions. In the large limit of Rd, 
we have confirmed that the fluid velocity converges on ui (77). Note 
that the cloud evaporates when Ud is positive, and vice versa. 

To estimate the front velocity, we need to know the relation¬ 
ship between Ud{R) and Rd- Here we evaluate / = —Rd/ud{R) 
from the results for d — 1 case. Noting that Cm in ea. <13> corre¬ 
sponds to the mass flux, Cm = —pu\, ui (77) can be written as 

Ul(R) = -CmXf^^ + °^^/p, ( 21 ) 

similarly, 

77i = -Mi(0) = c^X(0)^/<^+“Vp. (22) 

where the argument 0 represents values at r = 0. From those, 

/ = [X(0)/A'2]^/‘^+“^ . (23) 

Thus, the motion of the front is described by 
• d77 • .1 d — 1 

Rd = ^ = Ri-f--^X2. (24) 

dr p Rd 

FiglHshows the same as FigQbut for the front velocity Rd- At 
Rd ^ 10^, the approximate solutions well agree with the numerical 



10 100 1000 
Front Radius R 

Figure 1. Fluid velocity ltd (77) against radius 77 for (a) d = 3 (sphere) and 
(b) d = 2 (cylinder). The curves and symbols represent the approximate so¬ 
lution, eq. iini, and the numerical solution given by eqs. fT^ and fTTl . From 
the top to the bottom, p = 0.7, 0.8, 0.9,1,1.1,1.2 and 1.3, respectively. 
We set Z 2 = 2 and A = D = 1 for simplicity. 


ones. On the other hand. At Rd < 10^, the deviation of the approx¬ 
imate solutions from the numerical ones becomes large. This might 
suggest that the way of estimating / is too simple because we have 
found that the fluid velocity itself is well approximated by ea. <20t . 
Nevertheless, eQ. <24t well describes the critical radius, 77crit,d, at 
which Rd — 0, 


77 


crit,(i 


d- 1 

Cm 




d- 1^0 

pui 


(25) 


where T 2 = In fact, we have found that we can obtain 

a better fit when we replace / by while it makes the good 

estimation of 77crit,d worse. 

In the expression of eQ. <24> . an undesirable dependence on the 
temperature at the front remains, whose definition also has an ambi¬ 
guity. This should be removed by replacing X 2 and 77 in the curva¬ 
ture term by values at the maximum of X'X. To do this, however, 
the thickness of the front should be explicitly considered. This will 
be done in a subsequent work. Nevertheless, we have found that it 
provides a good fit to numerically obtained solutions when the front 
is defined at the unstable equilibrium in the case of the saturation 
pressure, p = 1. Apart from this, we would like to stress that the 
above approximate solution can be derived from values 77i and / 
in d = 1 case. 

Finally, we discuss the dependence of the front velocity on the 
front position. As mentioned in Section l. lShaviv & Regevi|il994h 
considered the frontal motion based on a model equation similar to 
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Figure 2. The same as Fig0 but for the front velocity, Rd- The curves 
and symbols represent the approximate solution, eQ. I24l . and the numer¬ 
ical solution given by eqs. d and From the top to the bottom, 

p = 1.3,1.2,1.1,1, 0.9, 0.8 and 0.7, respectively. 



Front Radius R 

Figure 3. The absolute value of the front velocity |i?| against R. To see its 
radial dependence, the vertical axis is logarithmic. The straight lines and 
symbols represent the front velocity provided by the approximate solution 
and the numerical solution, respectively. The upper and lower sets indicate 
d = 3 and 2, respectively. It is evidenUthatmjrresi^s^eproportional to 
R~^ different from the prediction by|Shaviv^Rege3|l994|). 


the TDGL equation. They found that R oc R~‘^ when pressure is 
nearly equal to the saturation pressure, and R oc when 

far from the saturation pressure. Because we have already shown 
in Figl^that R const, as R oo when far from the saturation 
pressure, we show the case of p = 1, in which we can see only the 
curvature effects. Figl^shows a log-log plot for |7?ci| against R for 
d = 2 and 3. It is evident that the front velocity is proportional to 
R~^ independent of d. Thus we conclude that the TDGL equation 
provided in the Lagrangian coordinate does not provide a correct 
model for thermally bistable fluids. 


4 GENERAL CURVED FRONTS 

Using the curvature term derived in the previous secti on, we discus s 
the dynamics of general curved fronts according to lBravl il994l . 
We define the x-axis normal to the unperturbed (straight) front, 
and take a unit vector normal to the curved front, g, with a di¬ 
rection from the CNM to the WNM. The situation under consid¬ 
eration is shown in Fig|4| schematically. Using this, we can write 
VX = {dgX)^g and = (dgX)^ -F (9gX)^ V • g. Substitut¬ 
ing these into equation <19L we obtain 

UdidgX)^ = [{d^^X)^ + {dgX)^V-g]j 

(26) 

a 

As done in the previous section, substituting corresponding terms 
in the right hand side with the one-dimensional equation and not¬ 
ing the existence of the direction cosine with respect to the normal 
direction of the unperturbed front, we obtain 

Vd = cose[Ri- fX 2 K/p\, (27) 

where Vd is the velocity of the front along the a;-axis, and K = X g 
is the mean curvature. When we take Rd as a curvature radius, we 
obtain K = {d— 1)/Rd. This result is applicable to general curved 
fronts. 

Let us consider a convex region of the CNM against the 
WNM, in which K > 0. When the CNM is evaporating (p < 1), 
the sign of Ri is negative. Because the sign of the second term 
is always negative, the CNM in the convex region evaporates faster 
than in concave regions. When the CNM is accreting cooling WNM 
(p > 1), the sign of Ri is positive. Thus the condensation in the 
convex region is slower then in concave regions. Consequently, the 
curvature term always smooths out the curved front. Note that this 
conclusion is valid only under the assumptions that flows are nor¬ 
mal to the front and that the structure of the front is the same as that 
in the case of plane-parallel geometry, in addition to the approxi¬ 
mation of isobaric evolution. 


5 CONCLUSIONS AND DISCUSSION 

We have investigated the dynamics of thermally bistable fluids from 
a pattern-theoretical point of view. To evaluate the curvature effects 
of the front connecting the WNM and CNM, at first, we focused on 
d-dimensional spherically symmetric CNM clouds. Using a way of 
approximations often used in the field of pattern formation theories, 
we derived an approximate solution describing the velocity of the 
front. Comparing with numerical solutions, we confirmed that they 
are in good agreement with each other, at least when the cloud size 
is larger than a few tens times the Field length, which corresponds 
to the thickness of the front. We have also found that our results 































CNM 





Figure 4. Schematic view of a curved front. The thick wavy curve rep¬ 
resents the front connecting the CNM {left) and WNM {right). The solid 
arrows denote the normal vector g, and the dashed arrow the velocity of the 
front along the horizontal axis . 


contradict those given bv ishaviv & Rese^ il994h . which assumed 
a model equation similar to the TDGL equation. We showed that 
the velocity of the front is proportional to the inverse of the radius 
in the case of the saturation pressure, and is constant in the case 
of much larger clouds and/or pressure far from the saturation pres¬ 
sure. Second, we discussed the dynamics of general curved fronts. 
Using the obtained approximate solution, we have found that the 
curvature effects smooth out curved fronts. 

In contrast to the latter concl usion, recent numerical ex- 

S ents of interstellar turbulence iKovama & Inutsu^ 120029. 

show that most fronts might be unstable even without 
strong shocks. This apparent contradiction might come from 
our simple assumptions that the structure of the front is in¬ 
dependent of geometry, that the fluid motion is normal to 
the front, and that the fluid evolves isobarically. Presumably 
some instability mechanisms are there similar to the Darrieus- 
Landau instability in propagating flame jUandau & Lifshit^l 19871 : 
llnoue. Inutsuka & Kovamal2005l). or the Mullins-Sek erka instabil- 
itv in crystal growth iMullins & S ekerkal 1 96.31 1 1 963) . This work is 
a first step toward full understanding of the dynamical behaviour 
of the ISM found in the numerical experiments. For fair compari¬ 
son, we must clarify at least how the above assumptions affect the 
dynamics of the ISM. We are planning to analyse the stability of 
the front considering these known mechanisms, as well as the va¬ 
lidity of the isobaric evolution. Together with those, hydrodynamic 
simulations should be performed to be compared with the results 
obtained in this paper. 
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